This report evaluates the simulation results.
Table showing the results between whether the transcript was set to be differentially expressed (DE) and if it overlaps (minimum 1 bp) any candidate DER.
## Overlaps DER
## DE status FALSE TRUE Sum
## FALSE 25 23 48
## TRUE 5 43 48
## Sum 30 66 96
The results are not the same using a minimum overlap of 13 bp between transcripts and candidate DERs with a FWER adjusted p-value < 0.05. Thus, we will use only the DERs with FWER adjusted p-value < 0.05.
## Overlaps DER (sig FWER)
## DE status FALSE TRUE Sum
## FALSE 37 11 48
## TRUE 6 42 48
## Sum 43 53 96
At a finer level, there is a difference in the number of exons per transcript overlapping all candidate DERs vs the DERs with FWER adjusted p-value < 0.05.
##
## 0 1 2 3 4 6 7 8 9 11 12 15 16 17 53 64 121 125
## 48 11 2 9 4 5 2 1 2 3 2 1 1 1 1 1 1 1
We can separate the transcripts by their experiment setup case. That is, whether its from a gene with:
Then compare against the results where
| Failed.DE | Failed.DER | Success.DE | Success.DER | |
|---|---|---|---|---|
| bothDE | 0 | 0 | 24 | 0 |
| noneDE | 0 | 0 | 0 | 24 |
| oneDE | 1 | 11 | 11 | 1 |
| singleDE | 5 | 0 | 7 | 0 |
| singleNotDE | 0 | 0 | 0 | 12 |
The 6 Failed.DE cases (false negatives) are mostly short single transcript genes (one exon only) where 5 were set to have low expression on one group, normal on the other two.
| tx_idx | tx_n | tx_i | gene_id | ucsckg_id | fasta_i | DE | group1 | group2 | group3 | width | readspertx | mean1 | mean2 | mean3 | case | result | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2 | 2 | 1 | 1 | 100126318 | uc021wmk.1 | 181 | TRUE | low | normal | normal | 78 | 31 | 16 | 31 | 31 | singleDE | Failed.DE |
| 11 | 13 | 1 | 1 | 100302149 | uc021wrh.1 | 796 | TRUE | normal | normal | low | 66 | 26 | 26 | 26 | 13 | singleDE | Failed.DE |
| 17 | 20 | 1 | 1 | 100422998 | uc021wnn.1 | 317 | TRUE | low | normal | normal | 86 | 34 | 17 | 34 | 34 | singleDE | Failed.DE |
| 22 | 27 | 1 | 1 | 100500833 | uc010gvn.2 | 347 | TRUE | normal | normal | high | 110 | 44 | 44 | 44 | 88 | singleDE | Failed.DE |
| 24 | 29 | 1 | 1 | 100500901 | uc021wny.1 | 423 | TRUE | normal | normal | low | 58 | 23 | 23 | 23 | 12 | singleDE | Failed.DE |
| 126 | 272 | 2 | 2 | 3761 | uc003avt.2 | 605 | TRUE | normal | normal | low | 2075 | 830 | 830 | 830 | 415 | oneDE | Failed.DE |
However, 2 similar cases with short transcripts were successfully detected. So it's likely that a lower F-stat cutoff would have picked up these false negative cases.
| tx_idx | tx_n | tx_i | gene_id | ucsckg_id | fasta_i | DE | group1 | group2 | group3 | width | readspertx | mean1 | mean2 | mean3 | case | result | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 10 | 12 | 1 | 1 | 100302118 | uc021wls.1 | 127 | TRUE | normal | high | normal | 78 | 31 | 31 | 62 | 31 | singleDE | Success.DE |
| 23 | 28 | 1 | 1 | 100500860 | uc021wlo.1 | 115 | TRUE | normal | low | normal | 88 | 35 | 35 | 18 | 35 | singleDE | Success.DE |
More info:
## IntegerList of length 6
## [["uc021wmk.1"]] 78
## [["uc021wrh.1"]] 66
## [["uc021wnn.1"]] 86
## [["uc010gvn.2"]] 110
## [["uc021wny.1"]] 58
## [["uc003avt.2"]] 220 1844
## IntegerList of length 2
## [["uc021wls.1"]] 78
## [["uc021wlo.1"]] 88
Coverage plots with F-statistics shown at the bottom for the false negative cases. One plot it shown for each exon that compose these transcripts.
## Matching regions to genes.
## nearestgene: loading bumphunter hg19 transcript database
## finding nearest transcripts...
## AnnotatingDone.
Out of the 11 Failed.DER transcripts (false positives), 11 of them are from the oneDE case. You could then argue that they are really not false positives. However, 0 and 0 transcripts are from the noneDE and singleNotDE cases respectively which would be the truly false positives.
| tx_idx | tx_n | tx_i | gene_id | ucsckg_id | fasta_i | DE | group1 | group2 | group3 | width | readspertx | mean1 | mean2 | mean3 | case | result | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 26 | 6 | 2 | 2 | 100130717 | uc011agh.3 | 17 | FALSE | normal | normal | normal | 650 | 260 | 260 | 260 | 260 | oneDE | Failed.DER |
| 37 | 56 | 2 | 1 | 10738 | uc010gwn.3 | 467 | FALSE | normal | normal | normal | 1488 | 595 | 595 | 595 | 595 | oneDE | Failed.DER |
| 47 | 80 | 2 | 1 | 128977 | uc002zpi.3 | 84 | FALSE | normal | normal | normal | 1558 | 623 | 623 | 623 | 623 | oneDE | Failed.DER |
| 49 | 85 | 2 | 1 | 129138 | uc003auc.3 | 576 | FALSE | normal | normal | normal | 2149 | 860 | 860 | 860 | 860 | oneDE | Failed.DER |
| 98 | 206 | 2 | 2 | 25817 | uc003bio.4 | 836 | FALSE | normal | normal | normal | 2652 | 1061 | 1061 | 1061 | 1061 | oneDE | Failed.DER |
| 121 | 266 | 2 | 1 | 339669 | uc003aqe.3 | 538 | FALSE | normal | normal | normal | 1049 | 420 | 420 | 420 | 420 | oneDE | Failed.DER |
| 132 | 283 | 2 | 2 | 3976 | uc011aks.2 | 379 | FALSE | normal | normal | normal | 3987 | 1595 | 1595 | 1595 | 1595 | oneDE | Failed.DER |
| 150 | 321 | 2 | 2 | 4689 | uc003apz.4 | 533 | FALSE | normal | normal | normal | 1646 | 658 | 658 | 658 | 658 | oneDE | Failed.DER |
| 175 | 400 | 2 | 1 | 6523 | uc003amc.3 | 457 | FALSE | normal | normal | normal | 5061 | 2024 | 2024 | 2024 | 2024 | oneDE | Failed.DER |
| 190 | 417 | 2 | 2 | 6948 | uc003air.2 | 407 | FALSE | normal | normal | normal | 2006 | 802 | 802 | 802 | 802 | oneDE | Failed.DER |
| 192 | 421 | 2 | 2 | 7122 | uc010grr.2 | 91 | FALSE | normal | normal | normal | 1720 | 688 | 688 | 688 | 688 | oneDE | Failed.DER |
Coverage plots with F-statistics shown at the bottom for the false positive cases. One plot it shown for each exon that compose these transcripts. For the 11 transcripts from the oneDE case, it can be seen how at least one plot contains a DER overlapping an exon set to be DE. .
Some complex situations where there are exons on both strands can be observed.
## Matching regions to genes.
## nearestgene: loading bumphunter hg19 transcript database
## finding nearest transcripts...
## AnnotatingDone.
In some simulations, we found what seemed to be false positive transcripts but turned out to overlap DERs in regions where there are exons on both the positive and negative strands and at least one of the exons was set to be DE.
## < table of extent 0 >
## IntegerList of length 0
## integer(0)
Skipped the following section due to the absence of such cases in this particular simulation
As it can be seen below, 0 apparent false positive transcripts from the noneDE case overlap (when strand is not taken into account) genes where at least one of two transcripts was set to be DE.
| tx_idx | tx_n | tx_i | gene_id | ucsckg_id | fasta_i | DE | group1 | group2 | group3 | width | readspertx | mean1 | mean2 | mean3 | case | result |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
For each gene, if at least one transcript is set to be DE then we consider the gene to be DE. Then, we check if the gene overlaps at least one DER.
## Overlaps DER
## DE status FALSE TRUE Sum
## FALSE 24 0 24
## TRUE 6 30 36
## Sum 30 30 60
The results from the simulation are promising as most transcripts were correctly classified as differentially expressed or not by derfinder.
The majority of the false negative cases involved short single transcript genes with one group having low expression relative to the other two. These cases could potentially be mitigated by lowering the F-statistic threshold used in the derfinder analysis.
In some simulations there are some apparent false positives which are due to transcripts on one strand set not to be DE overlapping transcripts from the other strand set to be DE. This situation could be solved with strand-specific RNA-seq data and running derfinder for each strand separately.
Minimum number of reads per transcript as well as per sample.
## Distribution of the minimum number of reads per transcript
summary(apply(readmat, 1, min))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.0 213.2 453.5 625.5 783.0 2747.0
## Distribution of the minimum number of reads per sample
summary(apply(readmat, 2, min))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.000 7.000 9.000 9.067 11.000 19.000
The minimum number of reads per transcript for a given sample is 1.
Next we can evaluate the simulation by classifying the exonic segments as whether they should be DE or not. Then, we can find if the DERs overlap such segments. We would expect that the DERs with a FWER adjusted p-value < 0.05 would only overlap segments that were set to be DE.
We can check how many segments each DER overlaps. Ideally they should all overlap at least one segment, but there are some cases where this could not happen (6 in this case). Mostly because of small mismatches between the transcripts and the actual mRNA used in the simulation.
## Do DERs overlap at least one segment?
table(countOverlaps(fullRegions, segs))
##
## 0 1 2 3 4
## 6 451 9 1 2
## Widths of DERs not overlapping any segment
width(fullRegions[is.na(fullRegions$overlap)])
## [1] 39 4 3 1 1 1
## Do these DERs have a significant FWER adjusted p-value?
table(fullRegions$significantFWER[is.na(fullRegions$overlap)])
##
## TRUE FALSE
## 0 6
However, the main result is whether the DERs overlap segments expected to be DE. Note that for this comparison, DERs are unstranded and could potentially overlap two segments from different strands where only one of them was set to be DE.
## FWER p-value < 0.05
## Overlaps a DE segment TRUE FALSE Sum
## FALSE 0 30 30
## TRUE 122 311 433
## Sum 122 341 463
Regardless of whether the DER p-value, we see that 6.48 percent of the DERs overlapping at least one segment, incorrectly overlap a segment set not to be DE.
None of the DERs with a FWER adjusted p-value < 0.05 overlaps a segment set not to be DE, which is what we would expect.
## [1] "2014-12-04 22:41:09 EST"
## user system elapsed
## 187.738 5.959 196.395
## setting value
## version R Under development (unstable) (2014-11-01 r66923)
## system x86_64, darwin10.8.0
## ui X11
## language (EN)
## collate en_US.UTF-8
## tz America/New_York
## package * version date source
## acepack 1.3.3.3 2013-05-03 CRAN (R 3.2.0)
## AnnotationDbi * 1.29.10 2014-11-22 Bioconductor
## base64enc 0.1.2 2014-06-26 CRAN (R 3.2.0)
## BatchJobs 1.4 2014-09-24 CRAN (R 3.2.0)
## BBmisc 1.7 2014-06-21 CRAN (R 3.2.0)
## Biobase * 2.27.0 2014-10-14 Bioconductor
## BiocGenerics * 0.13.2 2014-11-18 Bioconductor
## BiocParallel 1.1.9 2014-11-24 Bioconductor
## biomaRt 2.23.5 2014-11-22 Bioconductor
## Biostrings 2.35.5 2014-11-24 Bioconductor
## biovizBase 1.15.0 2014-10-14 Bioconductor
## bitops 1.0.6 2013-08-17 CRAN (R 3.2.0)
## brew 1.0.6 2011-04-13 CRAN (R 3.2.0)
## BSgenome 1.35.8 2014-11-21 Bioconductor
## bumphunter * 1.7.2 2014-11-19 Bioconductor
## checkmate 1.5.0 2014-10-19 CRAN (R 3.2.0)
## cluster 1.15.3 2014-09-04 CRAN (R 3.2.0)
## codetools 0.2.9 2014-08-21 CRAN (R 3.2.0)
## colorspace 1.2.4 2013-09-30 CRAN (R 3.2.0)
## DBI 0.3.1 2014-09-24 CRAN (R 3.2.0)
## derfinder * 1.1.14 2014-11-22 Github (lcolladotor/derfinder@24f9fbb)
## derfinderHelper * 1.1.5 2014-11-05 Bioconductor
## derfinderPlot * 1.1.5 2014-11-05 Bioconductor
## devtools * 1.6.1 2014-10-07 CRAN (R 3.2.0)
## dichromat 2.0.0 2013-01-24 CRAN (R 3.2.0)
## digest 0.6.4 2013-12-03 CRAN (R 3.2.0)
## doRNG 1.6 2014-03-07 CRAN (R 3.2.0)
## evaluate 0.5.5 2014-04-29 CRAN (R 3.2.0)
## fail 1.2 2013-09-19 CRAN (R 3.2.0)
## foreach * 1.4.2 2014-04-11 CRAN (R 3.2.0)
## foreign 0.8.61 2014-03-28 CRAN (R 3.2.0)
## formatR 1.0 2014-08-25 CRAN (R 3.2.0)
## Formula 1.1.2 2014-07-13 CRAN (R 3.2.0)
## GenomeInfoDb * 1.3.7 2014-11-15 Bioconductor
## GenomicAlignments 1.3.12 2014-11-27 Bioconductor
## GenomicFeatures * 1.19.6 2014-11-04 Bioconductor
## GenomicFiles 1.3.8 2014-11-12 Bioconductor
## GenomicRanges * 1.19.18 2014-12-01 Bioconductor
## GGally 0.4.8 2014-08-26 CRAN (R 3.2.0)
## ggbio 1.15.0 2014-10-14 Bioconductor
## ggplot2 1.0.0 2014-05-21 CRAN (R 3.2.0)
## graph 1.45.0 2014-10-14 Bioconductor
## gridExtra 0.9.1 2012-08-09 CRAN (R 3.2.0)
## gtable 0.1.2 2012-12-05 CRAN (R 3.2.0)
## Hmisc 3.14.5 2014-09-12 CRAN (R 3.2.0)
## htmltools 0.2.6 2014-09-08 CRAN (R 3.2.0)
## IRanges * 2.1.24 2014-12-01 Bioconductor
## iterators * 1.0.7 2014-04-11 CRAN (R 3.2.0)
## knitr * 1.7 2014-10-13 CRAN (R 3.2.0)
## knitrBootstrap 1.0.0 2014-11-03 Github (jimhester/knitrBootstrap@76c41f0)
## lattice 0.20.29 2014-04-04 CRAN (R 3.2.0)
## latticeExtra 0.6.26 2013-08-15 CRAN (R 3.2.0)
## locfit * 1.5.9.1 2013-04-20 CRAN (R 3.2.0)
## markdown 0.7.4 2014-08-24 CRAN (R 3.2.0)
## MASS 7.3.35 2014-09-30 CRAN (R 3.2.0)
## Matrix 1.1.4 2014-06-15 CRAN (R 3.2.0)
## matrixStats 0.10.3 2014-10-15 CRAN (R 3.2.0)
## mime 0.2 2014-09-26 CRAN (R 3.2.0)
## munsell 0.4.2 2013-07-11 CRAN (R 3.2.0)
## nnet 7.3.8 2014-03-28 CRAN (R 3.2.0)
## OrganismDbi 1.9.0 2014-10-14 Bioconductor
## pkgmaker 0.22 2014-05-14 CRAN (R 3.2.0)
## plyr 1.8.1 2014-02-26 CRAN (R 3.2.0)
## proto 0.3.10 2012-12-22 CRAN (R 3.2.0)
## qvalue 1.41.0 2014-10-14 Bioconductor
## R.methodsS3 1.6.1 2014-01-05 CRAN (R 3.2.0)
## RBGL 1.43.0 2014-10-14 Bioconductor
## RColorBrewer 1.0.5 2011-06-17 CRAN (R 3.2.0)
## Rcpp 0.11.3 2014-09-29 CRAN (R 3.2.0)
## RCurl 1.95.4.3 2014-07-29 CRAN (R 3.2.0)
## registry 0.2 2012-01-24 CRAN (R 3.2.0)
## reshape 0.8.5 2014-04-23 CRAN (R 3.2.0)
## reshape2 1.4 2014-04-23 CRAN (R 3.2.0)
## rmarkdown * 0.3.3 2014-09-17 CRAN (R 3.2.0)
## rngtools 1.2.4 2014-03-06 CRAN (R 3.2.0)
## rpart 4.1.8 2014-03-28 CRAN (R 3.2.0)
## Rsamtools 1.19.11 2014-11-26 Bioconductor
## RSQLite 1.0.0 2014-10-25 CRAN (R 3.2.0)
## rstudioapi 0.1 2014-03-27 CRAN (R 3.2.0)
## rtracklayer 1.27.6 2014-11-26 Bioconductor
## S4Vectors * 0.5.11 2014-11-24 Bioconductor
## scales 0.2.4 2014-04-22 CRAN (R 3.2.0)
## sendmailR 1.2.1 2014-09-21 CRAN (R 3.2.0)
## stringr 0.6.2 2012-12-06 CRAN (R 3.2.0)
## survival 2.37.7 2014-01-22 CRAN (R 3.2.0)
## TxDb.Hsapiens.UCSC.hg19.knownGene * 3.0.0 2014-09-26 Bioconductor
## VariantAnnotation 1.13.12 2014-11-20 Bioconductor
## XML 3.98.1.1 2013-06-20 CRAN (R 3.2.0)
## xtable 1.7.4 2014-09-12 CRAN (R 3.2.0)
## XVector * 0.7.3 2014-11-24 Bioconductor
## yaml 2.1.13 2014-06-12 CRAN (R 3.2.0)
## zlibbioc 1.13.0 2014-10-14 Bioconductor